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Abstract 

We have apphed the Finite Element Method to the seh-consistent elec- 
tronic structure calculations of molecules and solids for the first time. In 
this approach all the calculations are performed in " real space" and the use 
of non-uniform mesh is made possible, thus enabling us to deal with local- 
ized systems with ease. To illustrate the utility of this method, we perform 
an all-electron calculation of hydrogen molecule in a supercell with LDA 
approximation. Our method is also applicable to mesoscopic systems. 



I. Introduction 



Electronic-structure calculations of materials have a long history and var- 
ious methods have been developed for that purpose. One of the most popular 
methods for electronic-structure calculations of condensed matter systems today 
is the use of plane-wave basis sets and pseudo-potentials (PWPP) with local 
density approximation (LDA) [jl|. This method usually has enough accuracy for 
quantitative discussion, and reasonable efficiency to treat large systems. How- 
ever, for localized systems such as the first row elements or transition metals, the 
use of PWPP leads to a large number of plane-waves, thus making the calculation 
almost impossible. 

One way to avoid this difficulty is to use the ultrasoft pseudo-potentials 0. 
In this method, the constraint of norm-conservation is relaxed, and the number of 
plane-waves necessary for the calculation is considerably reduced. However, the 
deficit charge must be restored at every step, and the calculation becomes compli- 
cated. Furthermore, for systems with large vacuum regions such as molecules or 
slab model of surfaces in supercell geometry, even the ultrasoft pseudo-potential 
does not significantly reduce the number of plane-waves. So an efficient method 
to deal with such systems is desirable. 

II. Details of the Calculation 

Here we propose another way to solve above difficulties. In our method, all 
the calculations are performed in real space, and non-uniform meshes are used. 
First, the unit cell is divided into ortho-rhombic elements in real space, and the 
values of the wavefunctions on each vertex of the elements are taken as the basic 
variables. Then the wavefunctions are interpolated in each element as follows. 
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where hi,h2, denote the length of the three sides of the element and ipi, - ■ ■ ,ip8 
denote the value of the wavefunction on each vertex of the element (Fig.l). Note 
that u>s equals 1 at the s-th vertex and at the other vertices. 

In this way, the wavefunctions arc approximated by continuous piecewise poly- 
nomial functions. Here we put the values of the wavefunction of the i-th band on 
every vertex in the unit cell in some appropriate order and write it as ipi. 

The density corresponding to the wavefunction is approximated in the same 
way and represented by n. Here the values of the density is calculated from 
the square sum of the wavefunctions. In order to conserve the total charge, re- 
normalization of the density is performed. Though the present treatment of the 



density is not satisfactory, it causes no serious troubles such as instabilities, so 
we adopt it for the time being. 

In this approach, the total energy of the system is represented in the follow- 
ing way. First, we confine ourselves to the periodic systems to handle the long 
range Hartree potential efficiently. Furthermore, for simplicity we consider only 
real wavefunctions and bare coulomb potentials, though the extension to com- 
plex wavefunctions and separable norm-conserving pseudo-potentials is straight- 
forward. 

The total energy per unit cell is given by 

Etotal — Ekin + Epot + Ehartree + E^c + Eg^jald (10) 

in the LDA approximation. Rydberg units are used throughout the paper. 
The explicit form of the kinetic energy is 

Ekin^Y.! (11) 
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This integral can be calculated exactly in each element for the interpolated 
expression, and the result is the quadratic form of the values of wavefunctions at 
the vertices. So, with the aid of a constant matrix G, we can write 

Ek^n = Y.''Pi-G-^i- (12) 

i 

We now go on to the potential energy. After the divergences of the long range 
Coulomb potentials of nuclei and electrons are cancelled, the potential energy is 
expressed as 

Epot = -Stt E ^ ■ (^E Z^ exp (^ G ■ f,) j , (13) 

where n{G) is the reciprocal-space representation of the density and and fj 
denote the charge and the position of each nucleus in the unit cell respectively. 



Here we define a potential in real space as 

V (r-) = --^ E E ^ (5 ■ (f - rl)) ) . (14) 

cell i \g^q J 

In practice, the above sum is calculated with Ewald's method 0. Using this 
potential, the potential energy can be written as 

Epot = I V{r)n{r)df. (15) 

J cell 

This can again be calculated in each element, and the result is the linear form of 
the density. So we can write with the use of a constant vector 

Epot = 'v-n, (16) 

where each element of v is calculated as the sum of several integrals of the form 

lelement V (r) df. 

In the same way, the Hartree energy can be written as 

= 4./ ^(,>(r>/f. (17) 
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with 

^{r-) = T.^^^Pi^G-r). (18) 

Note that Lp (r) depends on the desity in contrast with V {•r). As a result, (f (r) 
must be calculated every time the density is changed. The direct use of the right 
hand side of (|18D is costly, because the fast Fourier transform (FFT) is necessary 
to obtain the values of n (G) for all G. So we calculate Lp (r) by minimizing 

/ M = £^ {I - (n (f) -n)-ip (r)} c/r (19) 

with the conjugate gradient method under the contraint, 

ip{f)df=0, (20) 

cell 



where n denote the average of the density. It can be easily shown that (f (r) 
given in Eq.(|TB|) minimizes the functional (|T^) under the same constraint. This 
technique for calculating Hartree potential is adopted in for a different reason. 
Here we approximate ip (r) in the same way as the density, and write it as . In 
this method, need not be so exact when the total energy is far from convergent. 
Accordingly, when the density is changed, only a few iterations are necessary to 
get the new Lp . Similar idea on this procedure is described in 0. 

Now we can calculate the Hartree energy and the result is the bilinear form 
of n and . Thus we can write with a constant matrix F , 

Ehartree = 4:7T*lf- F ■ n. (21) 

The rest of the total energy, E^c and Eewaid can be easily handled, and will 
not be given here. 

Now that the total energy is expressed as the function of {"i/^j} , we go on to 
minimize it with the conjugate gradient method [|l| under constraints 

J ilj,ir}^jir)dr = 5ij, (22) 

or in our notation, 

'tfi- F-tfj = Sij. (23) 

Here F is a sparse and symmetric matrix whose diagonal element is propor- 
tional to the volume of the element to which the vertex belongs. So if we apply 
the conjugate gradient method directly when we are working with non-uniform 
meshes, the conjugate gradient vector is destroyed in the orthogonalization pro- 
cedure, and the minimum of the total energy is never reached. In order to avoid 
this difficulty, we transform the wavefunctions as 

^[ = T-^i forVi. (24) 



Then the constraints become 



(25) 



with 




(26) 



We choose the matrix T that makes F' nearly the unit matrix. Another possible 
solution is to use the unconstrained minimization @ , in which no explicit orthog- 
onalization is required. Note that this transformation is different from the usual 
pre-conditioning, which is determined by the form of the total energy. 



To demonstrate the utihty of the present method, we calculated the equilib- 
rium structure of hydrogen molecule in a cubic supercell as illustrated in Fig.2. 
An all-electron calculation was performed with non-uniform meshes shown in 
Fig. 3. At the boundary between elements with different size, the wavefunctions 
are constrained to be continuous. The calculated density is shown in Fig. 4. The 
values of the calculated bond length and vibrational frequency are given in Table 
I and compared with those from experiments and other theories. The agreement 
is satisfactory. 



III. Example 



Bond Length (a.u.) Vibrational Frequency [cm ^) 



This work 
Other theory 
Experiment 



1.46 4424 
1.45 4277 
1.40 4400 



Table I. The other theory is from Ref |10|, in which LDA calculations are 
performed with Gaussian orbitals. The same form for e^c is adopted in these 
calculations [ill, |12 . 



IV. Discussion 

The advantages of this method are as follows, (i) Since all calculations are 
performed in real space, no FFT is necessary in contrast with PWPP. (ii) We can 
use non-uniform meshes. In the terminology of PWPP, we can effectively change 
the cutoff-energy locally. Accordingly, we can express localized orbitals without 
any serious difficulties. Furthermore, we can reduce the number of variables 
considerably when we calculate the electronic-structure of surfaces in a supercell 
configuration as illustrated in Fig.5. Similar idea is introduced in Ref [^, ^, in 
which plane-wave is used with adative Riemannian metric, (iii) In this method, 
the effect of the potential of the nucleus appears only in an integral form, so 
the divergence of the bare coulomb potential at the origin does not cause any 
difficulty, (iv) Orbitals localized in given regions of space [|| and their overlap 
matrix are naturally treated in this approach, (v) Our method is much easier to 
implement than PWPP. 

One clear limitation of this method is that it will be ineffective if we perform 
dynamical simulations with non-uniform meshes, which lead to the Pulay-like 
forces 0. 

V. Conclusion 

In conclusion, we have applied the Finite Element Method to the self- 
consistent calculation of electrons successfully for the first time. Since this method 
has several advantages over PWPP, it will be considerably favorable for the calcu- 
lations of some systems. Extensions of this method to higher order interpolation 
formulas and separable pseudo-potentials are under way. 
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Fig.l. An element and its vertices. 



Fig. 2. Hydrogen molecule in a cubic supercell (12 a.u. on a side). The cal- 
culations were performed with an eighth of the supercell. 

Fig. 3. In practice, we used the mesh twice as dense as this figure. The mesh is 
taken approximately logarithmic. 

Fig.4. The calculated density of hydrogen molecule. The singularity at the 
nucleus is well described. 



Fig. 5. In the vacuum region, only a small number of elements will suffice. 



